---
jupytext:
  cell_metadata_json: true
  formats: md:myst,ipynb
  text_representation:
    extension: .md
    format_name: myst
    format_version: 0.13
    jupytext_version: 1.16.4
kernelspec:
  display_name: Python 3 (ipykernel)
  language: python
  name: python3
---

```{code-cell}
:init_cell: true

import mmf_setup;mmf_setup.nbinit()
%pylab inline --no-import-all
```

# Travelling Waves

+++

We start with a brief review of Galilean Covariance.  For details see:

* [Galilean Covariance](http://swan.physics.wsu.edu/forbes/public/student_resources/galilean-covariance/)

+++

## Boosts

+++

Consider a solution to the GPE in 1D:

$$
  \I\hbar \dot{\psi}(x, t) = 
  -\frac{\hbar^2\psi''(x, t)}{2m}
  + g\abs{\psi(x, t)}^2\psi(x, t).
$$

Suppose we wish to study a solution that travels with velocity $v$.  I.e. something which should be static when considered as a function of $x_v = x - vt$.  As discussed in [Galilean Covariance], we may consider several transformations of the form:

$$
  \psi(x, t) = e^{\I\phi(x_v, t)}\psi_v(x_v, t)
$$

* Simple boost (notation $\psi_v(x_v, t)$):

  $$
    \hbar\phi(x_v, t) = 0, \\
    \I\hbar\dot{\psi}_v(x_v, t) = -\frac{\hbar^2\psi_v''(x_v, t)}{2m} 
    + \overbrace{\I\hbar v \psi_v'(x_v)}^{-v\op{p}\psi_v}
    + g\abs{\psi_v(x_v, t)}^2\psi_v(x_v, t).
  $$

* Full Galilean Transformation (notation with a capital $V=v$, $\psi_V(x_V, t)$):

  $$
    \hbar\phi(x_V, t) = mVx_V + \tfrac{m}{2}V^2t, \\
    \I\hbar\dot{\psi}_V(x_V, t) = -\frac{\hbar^2\psi_V''(x_V, t)}{2m} + g\abs{\psi_V(x_V, t)}^2\psi_V(x_v, t).
  $$

The full transformation makes the covariance of the original equation manifest - the form of the equation stays the same.  However, it does not simplify if the original equations are not Galilean covariant (e.g. if there is a modified dispersion due to a spin-orbit coupling.)

[Galilean Covariance]: http://swan.physics.wsu.edu/forbes/public/student_resources/galilean-covariance/

+++

### Arbitrary Dispersion

+++

The previous relationships can be expressed for arbitrary dispersion, which reduce to the previous equations when $E(p) = p^2/2m$:

$$
  \I\hbar \dot{\psi}(x, t) = 
  E(\op{p})\psi(x, t)
  + g\abs{\psi(x, t)}^2\psi(x, t),\\
  \psi(x, t) = e^{\I\phi(x_v, t)}\psi_v(x_v, t) = e^{\I\phi(x_V, t)}\psi_V(x_V, t).
$$

* Simple boost:

  $$
    \hbar\phi(x_v, t) = 0, \\
    \I\hbar\dot{\psi}_v(x_v, t) = \bigl(E(\op{p})-v\op{p}\bigr)\psi_v(x_v,t) 
    + g\abs{\psi_v(x_v, t)}^2\psi_v(x_v, t).
  $$

* Full Galilean Transformation:

  $$
    \hbar\phi(x_V, t) = mVx_V + \tfrac{m}{2}V^2t, \\
    \I\hbar\dot{\psi}_V(x_V, t) = \left(E(\op{p} + mV)-V\op{p} - \frac{mV^2}{2}\right)\psi_V(x_V, t) + g\abs{\psi_V(x_V, t)}^2\psi_V(x_V, t).
  $$

+++

## Madelung Transform: Quantum Hydrodynamics

+++

Implementing a Madelung transformation we obtain the hydrodynamic equivalent by equating real and imaginary parts (multiplying through by $-2\sqrt{\rho}$ and $\sqrt{\rho}/\rho$ in each case:

$$
  \psi(x, t) = \sqrt{\rho(x, t)}e^{\I\phi(x, t)}, \qquad u(x, t) = \tfrac{\hbar}{m} \phi'(x, t),\\
  \dot{\rho} + (\rho u)' = 0, \qquad
  \mu - \hbar\dot{\phi} = g\rho + \frac{m u^2}{2} 
  - \frac{\hbar^2}{2m}\frac{\sqrt{\rho}''}{\sqrt{\rho}}.
$$

Taking the gradient of the last equation, multiplying through by $\rho/m$ and using the first continuity equation to replace $\rho\dot{u} = (\rho u)_{,t} + u (\rho u)'$ we obtain:

$$
  (\rho u)_{,t} + \left(\rho u^2 + \frac{g\rho^2}{2m}\right)' = \frac{\rho}{m}\left( 
    \frac{\hbar^2}{2m}\frac{\sqrt{\rho}''}{\sqrt{\rho}}
  \right)'
  = \frac{\hbar^2}{4m}\bigl(\rho\, (\log\rho)''\bigr)'.
$$

In the first form, we can identify the roles of the force $F = -\nabla V$ which includes both the mean-field pressure and any external potential, and the effect of the quantum pressure:

$$
  (\rho u)_{,t} + (\rho u^2)' = \frac{\rho}{m}F,
  \qquad
  F = - \nabla\left(
    (g\rho + V_{\mathrm{ext}})
    - \frac{\hbar^2}{2m}\frac{\sqrt{\rho}''}{\sqrt{\rho}}
  \right).
$$

The left-hand-side is simply the covariant fluid derivative of the momentum $m \rho^{-1}\d \vect{p}/\d t$ so Newton's law is manifest. The second form is the hydrodynamic form used in Hoefer and El.

+++

### Moving Frame

+++

If we implement the full Galilean transformation, then the equations of motion are invariant, but if we implement the partial transform then the equations are slightly modified:


$$
  \dot{\rho}_v
  + \Bigl((u_v\overbrace{-v}^{\text{boost}})\rho_v\Bigr)' = 0, \qquad
  \mu \overbrace{+ \frac{mv^2}{2}}^{\text{boost}}
  - \hbar \dot{\phi}_v 
  = g \rho_v
  + m\frac{(u_v\overbrace{-v}^{\text{boost}})^2}{2}
  - \frac{\hbar^2}{2m}\frac{\sqrt{\rho_v}''}{\sqrt{\rho_v}}.
$$

In contrast, the full Galilean transformation also shifts $u_v \rightarrow u_v + v$ and the chemical potential $\mu \rightarrow \mu - mv^2/2$ so that the equations revert to their original form.

+++

Note: for pedagogical purposes, we derive the first equation explicitly.  The partial Galilean transformation just changes the variable $x \rightarrow x_v$:

$$
  x_v = x - vt, \qquad
  \rho_v(x_v, t) = \rho(x, t) = \rho(x_v + vt, t), \qquad
  u_v(x_v, t) = u(x, t) = u(x_v+vt, t),
$$

hence, the partials are related by

$$
  \dot{\rho}_v = \dot{\rho} + v\rho', \qquad
  \rho' = \rho_v', \qquad
  u_v' = u'.
$$

This follows from the following notation:

$$
  \dot{\rho} = \left.\pdiff{\rho(x, t)}{t}\right|_{x},\qquad
  \rho' = \left.\pdiff{\rho(x, t)}{x}\right|_{t}, \qquad
  \dot{\rho}_v = \left.\pdiff{\rho_v(x_v, t)}{t}\right|_{x_v}, \qquad
  \rho_v' = \left.\pdiff{\rho_v(x_v, t)}{x_v}\right|_{t}.
$$

Thus, the continuity equation is as advertised:

$$
  \overbrace{\dot{\rho}}^{\dot{\rho}_v - v \rho_v'} + \overbrace{(\rho u)'}^{(\rho_v u_v)'} =
  \dot{\rho}_v + \Bigl((u_v-v) \rho_v\Bigr)' = 0.
$$

+++

### Arbitrary Dispersion

+++

In the presence of arbitrary dispersion $E(p)$, the second quantum hydrodynamic equation (the Bernoulli equation) is significantly complicated outside of homogeneous matter (see [Khamehchi:2017]), however, the continuity equation remains the same.  Thus, in a simply boosted frame, we still have

$$
  \dot{\rho}_v + \Bigl((u_v-v)\rho_v\Bigr)' = 0.
$$


[Khamehchi:2017]: https://doi.org/10.1103/PhysRevLett.118.155301 'M.~A. Khamehchi, Khalid Hossain, M.~E. Mossman, Yongping Zhang, Thomas Busch, Michael McNeil Forbes, and Peter Engels, "Negative mass hydrodynamics in a spin-orbit--coupled Bose-Einstein condensate", prl 118(15), 155301 (2017) [1612.04055](http://arXiv.org/abs/1612.04055)'

+++

## Traveling Waves

+++

Consider a solution where the density profile moves but maintains its shape.  We can then determine the phase from the continuity equation:

$$
  \rho(x, t) = \rho_v(\overbrace{x - vt}^{x_v}), \qquad
  \dot{\rho} + (\rho u)' = \Bigr(\overbrace{\bigl(u_v(x_v) - v\bigr)\rho_v(x_v)}^{-C}\Bigr)' = 0,\\
  u(x, t) = u_v(x_v) = v - \frac{C}{\rho(x, t)} = v - \frac{C}{\rho_v(x_v)}.
$$

This is simply a re-derivation of the boosted continuity equation above.  The wave-function can thus be rendered time-independent with an appropriately chosen chemical potential $\mu$:

$$
  \psi(x, t) = \sqrt{\rho_v(x_v)}e^{\I[\phi(x_v) - \mu t/\hbar]} = \psi_v(x_v)e^{-\I\mu t/\hbar}, \qquad
  \phi'(x_v) = u_v(x_v) = v - \frac{C}{\rho_v(x_v)},\\
  \mu\psi_v(x_v) = \bigl(E(\op{p}) - v\op{p}\bigr)\psi_v(x_v)
    + g\abs{\psi_v(x_v)}^2\psi_v(x_v).
$$

In the case of the Galilean covariant GPE we can go further and use the full Galilean transformation:

$$
  \psi(x, t) = \sqrt{\rho_V(x_V)}e^{\I[\phi(x_V) - \mu t/\hbar]}e^{\I\bigl[mVx_V + \tfrac{m}{2}V^2t \bigr]/\hbar}
  = \psi_{V}(x_V)e^{\I\bigl[mVx_V + \bigl(\tfrac{m}{2}V^2 - \mu\bigr)t \bigr]/\hbar},\\
  \mu\psi_V(x_V) = -\frac{\hbar^2\psi_V''(x_V)}{2m} 
    + g\abs{\psi_V(x_V)}^2\psi_V(x_V).
$$

It might look like one can thus express traveling waves in a Galilean covariant theory such as the GPE as purely real stationary solutions in an appropriately moving frame, but this only describes a subset of possible solutions because the non-linear piece $g\abs{\psi_V}^2$ couples both real and imaginary parts.

This subset of solutions can be found by direct integration starting from the initial condition that $\psi_v(0)>0$ is the maximum of the wave with $\psi_V'(0) = 0$.  Since the wave must be convex at this point, we must have:

$$
  \frac{\hbar^2\psi_V''(0)}{2m \psi_V(0)} = g\psi_V^2(0) - \mu  \leq 0.
$$

We thus have a continuous family of solutions with $0 < \psi_V(0) < \sqrt{\mu/g}$.

+++

### Numerical Example

+++

Here we numerically solve this equation using the scipy IVP solver, then compare this with our exact solution as implemented in the `gpe.exact_solutions.TravellingWave` class.

```{code-cell}
from scipy.integrate import solve_ivp

class ODE:
    hbar = 1.0
    g = 1.0
    mu = 1.0
    m = 1.0
    
    @property
    def psi_max(self):
        if self.mu/self.g > 0:
            return np.sqrt(self.mu/self.g)
        else:
            return 1.0

    def f(self, x, y):
        psi, dpsi = y
        ddpsi = 2*self.m*(self.g*np.abs(psi)**2 - self.mu)*psi/self.hbar**2
        return (dpsi, ddpsi)

    def event(self, x, y):
        """Termination when dpsi=0 again"""
        if x <= 0:
            return -1
        psi, dpsi = y
        return dpsi

    event.terminal = True
    event.direction = -1
    
    def solve(self, fraction=0.9, max_step=0.1):
        x_max = 10.0
        psi0 = fraction*self.psi_max
        ivp = solve_ivp(self.f, (0, x_max), (psi0, 0), 
                        max_step=max_step, events=[self.event])
        return ivp


ode = ODE()
ivp = ode.solve(0.2)
x, psi = ivp.t, ivp.y[0]
plt.figure(figsize=(10,3))
ax1 = plt.subplot(121)
ax1.plot(x, psi)
ax1.set(xlabel="x", ylabel=r'$\psi$')
ax2 = plt.subplot(122)
ax2.plot(x, psi**2)
ax2.set(xlabel="x", ylabel=r'$n=\psi^2$')
L = ivp.t_events[0]
```

Here we construct the same solution from our exact solution.

```{code-cell}
# Note: there is some error in the period for longer waves.
from gpe.imports import *
from gpe.exact_solutions import TravellingWaves

n1 = (psi**2).max()
n0 = 1e-8  # Should be zero, but for numerical reasons we just make small.
s = TravellingWaves(n1=n1, n0=n0, g=ode.g, Lx=ivp.t[-1]/2.0)
s.plot()
plt.plot(ivp.t-s.Lx/2, ivp.y[0]**2)
print("psi_0/sqrt(mu/g) = {}".format(np.sqrt(n1/(s.get_mu().real/s.g))))
```

## Velocity of a Travelling Wave

+++

In the previous discussion, $v$ is the phase velocity of the traveling wave. This is the velocity with which the peaks of the wave travel.  *In contrast, the group velocity determines how fast a wave packet will move, which is a different concept.*  This velocity, however, is frame dependent, and in general, we would prefer to characterize this with respect to some specific frames.

* For waves with infinite wavelength on a non-zero background density (such as dark solitons, and some  bright solitons), a preferred frame is defined by the background at infinity which is taken to be at rest.
* For infinite-wavelength bright solitons with no background density, no preferred frame exists.  By boosting, one can realize these solutions with any velocity.
* For small amplitude oscillations (phonons), a preferred reference frame is defined by a stationary background density.

The question is how to define the preferred frame for large-amplitude, finite-wavelength modes.  One possibility is to define the frame with no net current:

$$
  \bar{p} = \frac{1}{L}\int_0^L\d{x}\; m \rho(x,t)u(x,t)
          = \frac{m}{L}\int_0^L\d{x_v} \left(v\rho_v(x_v) - C\right)
          = mv\bar{\rho} - m C = 0.
$$

With this definition, we have phase velocity:

$$
  v = \frac{C}{\bar{\rho}}.
$$

+++

# Hoefer and El

+++

We now consider the solution of Hoefer and El.  They set $\hbar=m=g=1$ so that we have the following units (in $d=3$ dimensions):

$$
  [\hbar] = MD^2/T, \qquad
  [g] = MD^{2+d}/T^2, \qquad
  [m] = M\\
  M = m = 1, \qquad
  mg/h^2 = D^{d-2}
  D = \left(\frac{mg}{\hbar^2}\right)^{1/(d-2)} = 1, \qquad
  T = \frac{m}{\hbar}\left(\frac{mg}{\hbar^2}\right)^{2/(d-2)} = 1,\\
  M = m = 1, \qquad
  D = \frac{mg}{\hbar^2} = 1, \qquad
  T = \frac{m^3g^2}{\hbar^5} = 1,\\
$$

$$
  \I\dot{\psi} = -\frac{1}{2}\pdiff[2]{\psi}{x} + \abs{\psi}^2\psi,\qquad
  \frac{\psi_V''}{2} = \left(\abs{\psi_V}^2 - \mu + \frac{V^2}{2}\right)\psi_V.
$$

Or, using the Madelung transform:

\begin{gather}
  \psi(x, t) = \sqrt{\rho(x, t)}e^{\I\phi(x, t)}, \qquad u(x, t) = \phi'(x, t),\\
  \dot{\rho} + (\rho u)' = 0, \qquad
  (\rho u)_{,t} + \left(\rho u^2 + \frac{\rho^2}{2}\right)' 
  = \frac{1}{4}\bigl(\rho\, (\log\rho)''\bigr)'. \tag{2.111}
\end{gather}

They express their solution in terms of four constants $r_1\leq r_2\leq r_3\leq r_4$ which we can related to some more physical parameters such as the amplitude $a$, minimum and maximum densities $\rho_{\min,\max}$, period $L$, and phase velocity $V$;

\begin{align}
  \DeclareMathOperator{\sn}{sn}
  \rho_\min &= (r_1-r_2-r_3+r_4)^2, \\
  V &= \frac{r_1+r_2+r_3+r_4}{2}, \tag{2.118}\\
  a = \rho_\max - \rho_\min &= (r_2-r_1)(r_4-r_3), \tag{2.117}\\
  \frac{1}{l^2} &= (r_4-r_2)(r_3-r_1),\tag{2.117}\\
  m = al^2 &= \frac{(r_2-r_1)(r_4-r_3)}{(r_4-r_2)(r_3-r_1)}, \tag{2.119}\\
  \rho(x, t) &= \rho_\min + a \sn^2\left(\frac{x_v}{l}, m\right), \tag{2.117}\\
  L &= 2lK(m), \tag{2.121}
\end{align}

+++

We disagree with their definition of $C$ which is:

$$
 C = \frac{1}{8}(-r_1-r_2+r_3+r_4)(-r_1+r_2-r_3+r_4)\overbrace{(r_1-r_2-r_3+r_4)}^{\sqrt{\rho_\min}}.
$$

Instead, we find

$$
  C = \frac{\sqrt{\rho_\min\rho_\max(1+\rho_\min l^2)}}{l}
    =\sqrt{(-r_1-r_2+r_3+r_4)^2 - 3a}
     \sqrt{(-r_1+r_2-r_3+r_4)^2 - \frac{3}{l^2}}
     (r_1 - r_2 - r_3 + r_4)
$$

+++

In our code, we would like to specify $\rho_\min$, $\rho_\max$, the period $L$ and the phase velocity $V$.  To do this, we must solve for the Jacobi parameter $m$:

$$
  L = 2lK(m) = 2\sqrt{\frac{m}{a}}K(m).
$$

```{code-cell}
import gpe.utils;reload(gpe.utils)
import gpe.bec;reload(gpe.bec)
from gpe.utils import evolve
import gpe.exact_solutions;reload(gpe.exact_solutions)
from gpe.exact_solutions import TravellingWaves

# v_p is the velocity of the constructed wave.
# v_x is the velocity of the frame.
s = TravellingWaves(Nx=64*2, Lx=10.0, n0=1.0, n1=1.1, v_p=0.5, v_x=0.5)
t_max = 10.0

for y in evolve(s, steps=500, dt_t_scale=0.5, t_max=t_max):
    plt.clf()
    y.plot()
```

The current is:

$$
  \vect{j} = \rho u = V\rho - C.
$$

In the moving frame:

$$
  (\rho_V u_V)' = 0, \qquad
  \mu + \frac{V^2}{2} = \rho_V + \frac{u_V^2}{2} 
  - \frac{1}{2}\frac{\sqrt{\rho_V}''}{\sqrt{\rho_V}}.
$$

Thus, if $\rho_V(x_v)$ is specified, then $u_V(x_V) = C/\rho_V(x_V)$ will ensure that the conservation equation is satisfied if:

$$
  \mu + \frac{V^2}{2} = \rho_V + \frac{C^2}{2\rho_V^2} 
  - \frac{1}{2}\frac{\sqrt{\rho_V}''}{\sqrt{\rho_V}}.
$$

The solution of Hoefer and El can be can be expressed as (Note: $k=\sqrt{m}$ for $\sn(z,k)$ in some CASs.  We use $\sn(z;m)$ and $K(m)$ here):

$$
  \psi(x, t) = \sqrt{\rho(x, t)}e^{\I\phi(x, t)}, \qquad
  u(x, t) = \phi'(x, t),\\
  m = al^2, \qquad
  \rho(x, t) = \rho_{\min} + a \sn^2\bigl(\frac{x_v}{l};m\bigr), \qquad
  x_v = x - Vt, \qquad
  u(x, t) = V - \frac{C}{\rho(x, t)}, \qquad
  \rho_{\max} = a + \rho_{\min}, \qquad
  C = \frac{\sqrt{\rho_{\min}\rho_{\max}(1+\rho_{\min}l^2)}}{l},\qquad
  L = 2lK(m), \qquad
  Q = \frac{2\pi}{L},\\
  \bar{\rho} = \frac{1}{L}\int_0^{L}\d{x}\;\rho(x, t)
             = \rho_{\min} + \frac{a}{m}
             = \left(1 - \frac{1}{m}\right)\rho_{\min} + \frac{1}{m}\rho_\max
             = \rho_{\min} + \frac{1}{l^2},\\
  V_p = \frac{C}{\bar{\rho}} 
      = \frac{\sqrt{\rho_{\min}\rho_{\max}(1+\rho_{\min}l^2)}}
             {l\rho_{\min} + \frac{1}{l}}
$$

+++

Note that, averaging over a complete period $L = 2lK(m)$ is:

$$
  \bar{\rho} = \frac{1}{L}\int_0^L \rho(x)\d{x}
  = \rho_\min + a\frac{1}{L}\int_0^L \sn^2(\tfrac{x}{l}, k)\d{x}
  = \rho_\min + \frac{a}{m}
  = \left(1-\frac{1}{m}\right)\rho_\min + \frac{1}{m}\rho_\max.
$$

+++

$$
  V_p = \frac{C}{\bar{\rho}} 
      = \frac{\sqrt{1+\rho l^2}}{l + \frac{1}{l\rho}}
$$

```{code-cell}
import gpe.utils;reload(gpe.utils)
import gpe.bec;reload(gpe.bec)
from gpe.utils import evolve
import gpe.exact_solutions;reload(gpe.exact_solutions)
from gpe.exact_solutions import BrightSoliton

s = BrightSoliton(Nx=64, Lx=10.0, 
                  v=0.0, # v=10.0
                  v_x=10.0)
t_max = s.basis.Lxyz[0]/abs(s.v_x)

for y in evolve(s, steps=500, t_max=t_max):
    plt.clf()
    y.plot()
```

## Reformulation

+++

Here we present a reformulation of the coefficients into more physically meaningful parameters.  We start with the 

$$
  x_v = x - vt, \qquad
  u(x_v) = v - \frac{C}{\rho(x_v)}.
$$

This form immediately satisfies the continuity equation as long as $\rho(x_v) = \rho(x-vt)$:

$$
  \rho_{,t} + (u\rho)_{,x} = -v\rho'(x_v) + (v\rho(x_v) - C)' = 0
$$

Next we have:

$$
  \rho(x_v) = \rho_\min + a\sn^2\left(\frac{x_v}{l}; m\right).
$$

```{code-cell}
%pylab inline --no-import-all
from gpe.exact_solutions import sn, K
m = 0.999999999
L = 2*K(m)
x = np.linspace(-L,L,100)
plt.plot(x, sn(x, m)**2)
```

```{code-cell}
import numpy as np
np.random.seed(1)
r = np.random.random(5)
r[4] = -r[1]-r[2]-r[3]
C = 1./8*(-r[1]-r[2]+r[3]+r[4])*(-r[1]+r[2]-r[3]+r[4])*(r[1]-r[2]-r[3]+r[4])
V = (r[1]+r[2]+r[3]+r[4])/2
m = (r[2]-r[1])*(r[4]-r[3])/(r[4]-r[2])/(r[3]-r[1])
k = np.sqrt(m)
a = (r[4]-r[3])*(r[2]-r[1])
l = 1./np.sqrt((r[4]-r[2])*(r[3]-r[1]))
rho_min = (r[1]-r[2]-r[3]+r[4])**2
rho_max = rho_min + a
V, np.sqrt(rho_min*rho_max*(1+rho_min*l**2))/l, C
```

## Phonon Limit

In the limit of small amplitude $a$ we have the following reduction:

$$
  a\rightarrow 0, \qquad
  \bar{\rho} \rightarrow \rho, \qquad
  C \rightarrow \frac{\rho}{l}\sqrt{1+\rho l^2}}
\epsilon = \rho_\max - \rho_\min = \epsil we obtain

In the phonon limit we should have

$$
  \omega(q) = \pm\sqrt{q^2(\rho_0 + q^2/4)}, \qquad
  V_p = \frac{\omega(q)}{q} = \pm\sqrt{\rho_0 + q^2/4}, \qquad
  V_g = \omega'(q) = \frac{q^2+2\rho_0}{\pm\sqrt{q^2 + 4\rho_0}}
$$

+++

Thus, one can specify the following four parameters (and an overall arbitrary initial phase):

* $\rho_\min$ and $\rho_\max$: These determine the amplitude $a = \rho_\max - \rho_\min \geq 0$ of the wave.
* $L$: The period of the wave can be chosen implicitly by setting $l$ which will determine the period through $m=al^2$ and $L = 2l K(m)$.
* $V$: The speed of the wave may also be arbitrarily chosen.  This may seem strange since the speed of the waves should follow from the dispersion relationship, however, this is simply a reflection of the fact that we can look at the system in an arbitrarily moving frame.
* The rest frame is defined as the frame in which there is no net momentum flow.  Thus, the average current should be zero:

  $$
    0 = \bar{\vect{j}} = \frac{1}{L}\int_0^L \vect{j}(x, t)\d{x}
    = \frac{1}{L}\int_0^L \rho(x, t)u(x, t)\d{x}
    = \frac{1}{L}\int_0^L \Bigl(V\rho(x, t) - C\Bigr)\d{x}\\
    V_p = \frac{CL}{\int_0^L \rho(x, t)\d{x}} = \frac{C}{\bar{\rho}}
  $$
  
  where $\bar{\rho}$ is the average density across the cell.  Note: this depiction will fail for bright solitons on the vacuum since $\bar{\rho} = 0$ but in this case, the notion of the fixed background has no meaning.

+++

## Soliton Limit
In the long wavelength limit $m \rightarrow 1$:

$$
  K(m) \rightarrow \ln \frac{4}{\sqrt{1-m}}, \qquad
  l = \frac{1}{\sqrt{a}}, \qquad
  C = \rho_\max\sqrt{\rho_{\min}},\qquad
  V_p = \frac{C}{\rho_\max} = \sqrt{\rho_{\min}}.
$$

Thus, the speed of the traveling waves is set by the *minimum* density in the long wavelength limit.  This is consistent, for example, with the speed of a grey soliton which becomes stationary when it becomes dark (i.e. when the minimum density in the soliton becomes zero.)

+++

$$
  \hbar\omega(q) = \pm \sqrt{\frac{\hbar^2q^2}{2m}\left(\frac{\hbar^2q^2}{2m} + 2\mu\right)}\\
  \omega(q) = \pm \abs{q}\sqrt{\rho_0 + \frac{q^2}{4}}\\
$$

+++

$$ L = 2\sqrt{\frac{m}{a}}K(m) $$

+++

$$
  \partial_t[e^{\I\phi}f(x-vt)] = \I\dot{\phi}e^{\I\phi}f(x-vt) -ve^{\I\phi}f'(x-vt)\\
  \partial_t[\psi] = \I\partial_t\phi \psi - v\psi'(x-vt)
$$

```{code-cell}
from gpe.imports import *
import gpe.exact_solutions; reload(gpe.exact_solutions)
from gpe.exact_solutions import TravellingWaves

s = TravellingWaves(v_p=10.0, v_x=10.0, n0=0.1, Nx=32*4, Lx=5)
s.plot()
t_max = s.basis.Lxyz[0]/abs(s.v_p)

for y in evolve(s, steps=400, t_max=t_max):
    plt.clf()
    y.plot()
```

```{code-cell}
from gpe.imports import *
#from mmfutils.math.special import ellipkinv
from scipy.special import ellipj, ellipk
import scipy.integrate
import scipy as sp
import gpe.bec

def sn(u, m):
    return ellipj(u, m)[0]

def K(m):
    return ellipk(m)


def get_QV(m=0.1, a=0.1, rho_min=1.0, N=256):
    rho_max = rho_min + a
    #m = a*l**2
    l = np.sqrt(m/a)
    L = 2*l*K(m)
    dx = L/N
    x = np.arange(N)*dx
    rho = rho_min + a*sn(x/l, m)**2
    C = np.sqrt(rho_min*rho_max*(1+rho_min*l**2))/l
    V = np.trapz(C/rho, x)/L
    Q = 2*np.pi/L
    return Q, V

rho_min = 1.0
c = np.sqrt(rho_min)
ms = 1 - 10**np.linspace(-10, -0.001, 1000)

plt.figure(figsize=(20,5))
ax1 = plt.subplot(121)
ax2 = plt.subplot(122)
for a in [0.01, 0.2, 1.0]:
    Q, Vp = np.array([get_QV(m, a=a, rho_min=rho_min) for m in ms]).T
    l, = ax1.plot(Q**2, Vp/c, '-+', label=a)
    Q_ = (Q[1:]+Q[:-1])/2
    Vp_ = (Vp[1:]+Vp[:-1])/2
    Vg_ = Q_*np.diff(Vp)/np.diff(Q) +Vp_ 
    ax1.plot(Q_**2, Vg_/c, '--', c=l.get_c())
    
    ax2.plot(Q, Q*Vp, '-+', label=a)

q_max = 2.0
q = np.sqrt(np.linspace(0, q_max**2, 100))
Vp = np.sqrt(rho_min + q**2/4)
w = Vp*q
Vg = (q**2 +2*rho_min)/np.sqrt(q**2 +4*rho_min)
l, = ax1.plot(q**2, Vp/c, label='phonon')
ax1.plot(q**2, Vg/c, '--', c=l.get_c())
ax1.set_xlim(0, q_max**2)
ax1.set_ylim(0.99, 1.5)
plt.sca(ax1)
plt.axis([0,1.5,1,1.5])
ax1.set_xlabel("$q^2$")
ax1.legend()

plt.sca(ax2)
ax2.plot(q, w, label='phonon')
ax2.set_xlabel("$q$")
ax2.set_ylabel("$\omega(q)$")
ax2.set_xlim(0, 1.25)
ax2.set_ylim(0, 2)
#plt.axis([0,0.5,0,0.5])
ax2.legend()


get_QV(m=0.9999999999999999, a=1.0)
```

```{code-cell}
1.15*2.0
```

```{code-cell}
from gpe.soc import u
L = 2.5/3.0*u.micron
g = 0.00077
m = 87*u.u
hbar = u.hbar
u_length = g/hbar**2*m
Q = 2*np.pi/(L/u_length)
rho_max = 10530.0/u.micron**3*u_length**3
rho_min = 452.0/u.micron**3*u_length**3
a = rho_max - rho_min
get_QV(m=0.999999992, a=a)
```

```{code-cell}
:init_cell: true

from gpe.imports import *
import scipy.optimize
from scipy.special import ellipj, ellipk
import scipy.integrate
import scipy as sp
import gpe.bec

def sn(u, m):
    return ellipj(u, m)[0]

def cn(u, m):
    return ellipj(u, m)[0]

def K(m):
    return ellipk(m)

def get_m(a, L):
    def f(m):
        return 2*np.sqrt(m/a)*K(m) - L

    m0 = m1 = 0
    f0 = f1 = f(m0)
    while f1 < 0:
        m1 = (m1 + 1)/2.0
        f1 = f(m1)
    return sp.optimize.brentq(f, m0, m1)
```

```{code-cell}
class State(gpe.bec.State):
    def __init__(self, Nx=256, cells_x=10, 
                 l=1.1757305, b=1.0, a=0.1, V=0.0):
        m = a*l**2
        L = 2*l*K(m) * cells_x
        mu = ((a+3*b)*l**2+1)/2/l**2
        C = np.sqrt((a+b)*b*(1+b*l**2)/l**2)
        dx = L/Nx
        gpe.bec.State.__init__(self, Nxyz=(Nx,), Lxyz=(L,), 
                               m=1.0, g=1.0, hbar=1.0)
        
        x = self.basis.xyz[0]
        rho = b + a*sn(x/l, m)**2
        u = V - C/rho
        psi = np.sqrt(rho)*np.exp(1j*sp.integrate.cumtrapz(u, x, initial=0))
        self.data[...] = psi
        # Hpsi = np.fft.ifft(k**2*np.fft.fft(psi))/2 + (abs(psi)**2-mu)*psi
        #plt.plot(x, abs(Hpsi))
        
    def get_Vext(self):
        return 0.0
        
s = State(Nx=256)

e = EvolverABM(s, dt=0.5*s.t_scale)
with NoInterrupt() as interrupted:
    while not interrupted:
        e.evolve(100)
        e.y.plot()
        display(plt.gcf())
        plt.close('all')
        clear_output(wait=True)
```

```{code-cell}
import scipy.integrate
import scipy as sp

l = 1.1858
l = 1.1757305
b = 1.0
a = 0.1
V = 0.0
m = a*l**2
L = 2*l*K(m)
mu = ((a+3*b)*l**2+1)/2/l**2
C = np.sqrt((a+b)*b*(1+b*l**2)/l**2)
N = 256
dx = 10*L/N
x = np.arange(N)*dx - L/2
k = 2*np.pi * np.fft.fftfreq(N, dx)
rho = b + a*sn(x/l, m)**2
u = V - C/rho
#plt.plot(x, rho)
psi = np.sqrt(rho)*np.exp(1j*sp.integrate.cumtrapz(u, x, initial=0))
Hpsi = np.fft.ifft(k**2*np.fft.fft(psi))/2 + (abs(psi)**2-mu)*psi
plt.plot(x, abs(Hpsi))
```

```{code-cell}
plt.plot(x, rho*u**2 + rho**2/2)
plt.plot(x[1:-1], 2.36+rho[1:-1]*np.diff(np.diff(np.log(rho)))/dx**2/4)
```

```{code-cell}
from gpe.imports import *
from gpe.bec import State
s = 
```

####

+++

## Phonons

+++

Phonons appear in the small amplitude limit $a\rightarrow 0$ in which case:

$$
  \sn(z;0) = \sin(z),\qquad
  K(m=0) = \frac{\pi}{2}, \qquad
  L = \pi l, \qquad
  Q = \frac{2}{l}\\
  \rho(x, t) = \rho_{\min} + a \sin^2\frac{Q x_v}{2}
             = \rho_{0} - \frac{a}{2}\cos(Q x_v), \qquad
$$

+++

$$
  C = \rho_0\sqrt{l^{-2}+\rho_{0}}\\
  u = V - \sqrt{\frac{Q^2}{4}+\rho_{0}}
$$

+++

The phase velocity is $V = \omega(Q)/Q$

+++

and present the density as:

$$
  \rho = \frac{1}{4}(r_4-r_3-r_2+r_1)^2 
    + u^2\sn^2\bigl(\frac{x_v}{l};m\bigr), \qquad
  x_v = x - vt\\
  v = \frac{n_1+n_2+n_3+n_4}{2}, \qquad
  a = u^2 = (r_4-r_3)(r_2-r_1), \\
  m = \frac{(r_2-r_1)(r_4-r_3)}{(r_4-r_2)(r_3-r_1)} 
    = u^2l^2, \qquad
  l = \frac{1}{\sqrt{(r_4-r_2)(r_3-r_1)}},\\
  c^2 = .
$$

The soliton limit is $m=1$ (i.e. $r_2=r_3$) where $\sn(z, 1) = \tanh(z)$ where the solution has the form:

$$
  \psi = \left[\I v + u \tanh\bigl(\frac{x_v}{l}\bigr)\right]e^{-\I (u^2+v^2) t}, \qquad
  u = 1/l,\\
  \rho = v^2 + u^2 \tanh^2\bigl(\frac{x_v}{l}\bigr)
       = v^2 + u^2 - u^2\sech^2\bigl(\frac{x_v}{l}\bigr)\\
  v^2 = \frac{1}{4}(r_4-2r_2+r_1)^2\\
  v^2 + u^2 = c^2 = \bar{\rho} = \frac{(r_4-r_1)^2}{4}.
$$

+++

$$
  (r_4-r_2)(r_2-r_1) + \frac{(r_4 - 2r_2 + r_1)^2}{4} = \frac{(r_1-r_4)^2}{4}
$$

+++

A traveling wave with velocity $v$ has the following form:

$$
  \psi(x, t) = \sqrt{\rho}(x - vt) e^{\I\phi(x-vt) - \I\omega t}, \\
  \I\hbar \dot{\psi} = -\frac{\hbar^2\psi''}{2m} + (g\rho - \mu) \psi.\\
  -\I\hbar v\psi' + \hbar\omega\psi = -\frac{\hbar^2\psi''}{2m} + (g\rho - \mu) \psi.  
$$

By appropriately adjusting $\mu$ we may set $\omega = 0$, hence we may drop the factor $e^{-\I\omega t}$ from the wavefunction:

$$
  \psi(x - vt) = \sqrt{\rho}(x - vt) e^{\I\phi(x-vt)}.
$$

The last equation above now represents a stationary solution in a moving frame with coordinate $y = x-vt$:

$$
  0 = -\frac{\hbar^2\psi''(y)}{2m} + \I\hbar v\psi'(y) + [g\rho(y)- \mu]\psi(y)
             = \left[\frac{\op{p}^2}{2m} -\op{p} v + g\rho(y) - \mu\right]\psi(y)
             = \left[\frac{(\op{p} - p_v)^2}{2m} - \frac{p_v^2}{2m} + g\rho(y) - \mu\right]\psi(y)
$$

where $p_v = mv$.  Finally, we may recast this in terms of the original GPE by transforming $\psi(y)$ as follows:

$$
  \psi(y) = e^{\I p_v x/\hbar}\tilde{\psi}(y)\\
  (\op{p}-p_v)\psi(y) = e^{-\I p_v x/\hbar}\op{p}\tilde{\psi}(y)\\
  \left[\frac{\op{p}^2}{2m} + g\rho(y) - \left(\mu + \frac{p_v^2}{2m}\right)\right]\tilde{\psi}(y)
$$

If the original function $\psi(y+L) = \psi(y)$ was periodic, then $\tilde{\psi}(y+L) = e^{\I p_v L/\hbar}\tilde{\psi}(y)$.  In other words, twisted boundary conditions should be applied with a twist $\theta = p_vL/\hbar$.

+++

# GPE

+++

$$
  \psi(x-vt) = e^{\I\phi(x-vt)}f(x-vt), \qquad f(x) = \sqrt{\rho(x)}\\
  \psi' = \I\phi'\psi + e^{\I\phi}f', \qquad
  \psi'' = \I\phi''\psi - (\phi')^2\psi +2\I\phi'\psi'
  + e^{\I\phi}\I\phi'f' + e^{\I\phi}f'',\\
  \I\dot{\psi} = -v\I\psi' = v\phi'\psi - v\I e^{\I\phi}f'
$$
$$
  2v\phi'\psi - 2v\I e^{\I\phi}f' = 
  \I\phi''\psi - (\phi')^2\psi +2\I\phi'\psi' + e^{\I\phi}\I\phi'f' + e^{\I\phi}f'' + 2e^{\I\phi}f^3\\
  [3\phi'+2v]f' = -\phi''f\\
  [3\phi' + 2v]\phi'f = f'' + 2f^3\\
  -\left(\frac{(\phi')^2}{2}\right)' 
  = \frac{f'f''}{f^2} + (f^2)'
  = \frac{f'f''}{f^2} + \rho'\\
$$

$$
  f^2 = \rho, \qquad
  2ff' = \rho', \qquad
  2f'f' + 2ff'' = \rho'', \qquad
$$

+++

Here we consider the problem of finding traveling wave solutions in a BEC.  From a numerical perspective, it is highly beneficial if the problem can be stated in terms of a well-defined minimization problem.  To start, we consider the available analytic solution for the conventional GPE.  These solutions are presented in [El:2016]:

$$
  \DeclareMathOperator{\sn}{sn}
  \I\dot{\psi} = \frac{-\psi''}{2} + \abs{\psi}^2\psi\\
  \rho = \abs{\psi}^2 = \frac{(r_4-r_3-r_2+r_1)^2}{4} 
  + (r_4-r_3)(r_2-r_1)\sn^2(\sqrt{(r_4-r_2)(r_3-r_1)}\xi; m)\\
  C = \frac{(-r_1-r_2+r_3+r_4)(-r_1+r_2-r_3+r_4)(r_1-r_2-r_3+r_4)}{8}\\
  u = v - \frac{C}{\rho}, \qquad
  \xi = x-vt - \xi_0, \qquad
  v = \frac{1}{2}\sum r_i, \\
  a = (r_4-r_3)(r_2-r_1), \qquad
  r_1 \leq r_2 \leq r_3 \leq r_4, \qquad
  m = \frac{(r_2-r_1)(r_4-r_3)}{(r_4-r_2)(r_3-r_1)},\\
  L = \frac{1}{2}\oint \frac{\d\lambda}
    {\sqrt{(\lambda - r_1)(\lambda - r_2)(\lambda - r_3)(\lambda - r_4)}}
    = \frac{2K(m)}{\sqrt{(r_4-r_2)(r_3-r_1)}}.
$$

The physical interpretations are:

* $a$: Amplitude
* $v$: Phase velocity
* $u$:

+++

$$
  \DeclareMathOperator{\sn}{sn}
  \I\dot{\psi} = \frac{-\psi''}{2} + \abs{\psi}^2\psi\\
  \rho = \abs{\psi}^2 = \alpha + a\sn^2(z; m)\\
  \alpha = \frac{(r_4-r_3-r_2+r_1)^2}{4}, \\
  D = \sqrt{(r_4-r_2)(r_3-r_1)}, \qquad
  z = D\xi\\
  C = \frac{(-r_1-r_2+r_3+r_4)(-r_1+r_2-r_3+r_4)(r_1-r_2-r_3+r_4)}{8}\\
  u = v - \frac{C}{\rho}, \qquad
  \xi = x-vt - \xi_0, \qquad
  v = \frac{1}{2}\sum r_i, \\
  a = (r_4-r_3)(r_2-r_1), \qquad
  r_1 \leq r_2 \leq r_3 \leq r_4, \qquad
  m = \frac{(r_2-r_1)(r_4-r_3)}{(r_4-r_2)(r_3-r_1)},\\
  L = \frac{1}{2}\oint \frac{\d\lambda}
    {\sqrt{(\lambda - r_1)(\lambda - r_2)(\lambda - r_3)(\lambda - r_4)}}
    = \frac{2K(m)}{\sqrt{(r_4-r_2)(r_3-r_1)}}.
$$

+++

$$
  \DeclareMathOperator{\cn}{cn}
  \DeclareMathOperator{\dn}{dn}
  \sn'(z) = \cn(z)\dn(z), \qquad
  \cn'(z) = -\sn(z)\dn(z), \qquad
  \dn'(z) = -k^2\sn(z)\cn(z)\\
  \sn''(z) = -\sn(z)[\dn^2(z) +k^2\cn^2(z)], \qquad
  \cn'(z) = -\sn(z)\dn(z), \qquad
  \dn'(z) = -k^2\sn(z)\cn(z)\\
  \cn^2 + \sn^2 = \dn^2 + k^2\sn^2 = 1.
$$

A special limit is when $m=1$:

$$
  \sn(z;m=1) = \tanh(z), \qquad
  \cn(z;m=1) = \dn(z;m=1) = \sech(z).
$$

+++

$$
  \psi = \sqrt{1+\alpha\sn^2(z)}\\
  \psi' = \frac{\alpha\sn(z)\sn'(z)}{\psi}\\
  \psi'' = \frac{\alpha}{\psi}\left(
  \sn'(z)\sn'(z) + \sn(z)\sn''(z)
  -\frac{\alpha\sn^2(z)[\sn'(z)]^2}{\psi^2}\right)\\
$$

+++

$$
  \psi' = \frac{\alpha\sn(z)\cn(z)\dn(z)}{\psi}\\
  \psi'' = \frac{\alpha}{\psi}\left(
  \cn^2(z)\dn^2(z) - \sn^2(z)\dn^2(z) - k^2\sn^2(z)\cn^2(z)
  -\frac{\alpha\sn^2(z)\cn^2(z)\dn^2(z)}{\psi^2}\right)\\
$$

+++

First we test these.  To match units we set $\hbar = m = g = 1$.

```{code-cell}
from gpe.imports import *
from scipy.special import ellipj, ellipk
import gpe.bec

def sn(u, m):
    return ellipj(u, m)[0]

def cn(u, m):
    return ellipj(u, m)[0]

def K(m):
    return ellipk(m)

class State(gpe.bec.State):
    def __init__(self, Nx=32, rs=[0.1, 0.2, 0.3, 0.4], xi0=0):
        r1, r2, r3, r4 = rs = sorted(rs)
        k = np.sqrt((r4-r2)*(r3-r1))
        m = (r2-r1)*(r4-r3)/(r4-r2)/(r3-r1)
        v = sum(rs)/2.0
        L = 2.*K(m)/k
        gpe.bec.State.__init__(self, Nxyz=(Nx,), Lxyz=(L,), 
                               m=1.0, g=1.0, hbar=1.0)
        x = self.xyz[0]
        self.T = abs(L/v)
        t = 0
        xi = x - v*t - xi0
        a = (r4-r3)*(r2-r1)
        #a = m**2*k**2
        rho_0 = (r4-r3-r2+r1)**2/4.0
        rho_0 = (r4-r3+r2-r1)**2/4.0
        rho = rho_0  + a*sn(k*xi, m)**2
        self.a_ = self.a = a
        self.m_ = m
        self.k_ = self.k = k
        self.rho_0_ = self.rho_0 = rho_0

        k_ = self.m*v/self.hbar
        self[...] = np.exp(1j*k_*x) * np.sqrt(rho)
        
    def get_Vext(self):
        return 0.0
        
s = State(Nx=256, rs=[-2.0, -1.0, 1.0, 2.0])
```

```{code-cell}
n = s.get_density()
x = s.xyz[0]
x_ = x[1:-1]
psi_ = s[1:-1]
ddpsi_ = np.diff(np.diff(s[...]))/np.diff(x)[1:]**2
plt.plot(x_, ddpsi_/psi_/2 + abs(psi_)**2)
```

Check the soliton limit (2.122)

```{code-cell}

```

```{code-cell}
s = State(Nx=128, rs=[-4.000001, 1.0, 1.000001, 2.0])
x = s.xyz[0]
s.plot()
rho_min = abs(s[...]).min()**2
rho_max = abs(s[...]).max()**2
plt.axhspan(rho_min, rho_min + s.a, fc='y', alpha=0.5)
#plt.plot(x, rho_max-s.a/np.cosh(np.sqrt(s.a)*x)**2, '+:')
```

```{code-cell}
rho_0 = s.rho_0
s[...] = np.sqrt(rho_0)*np.tanh(np.sqrt(rho_0)*x)
s.plot()
```

```{code-cell}
s.cooling_phase = 1
e = EvolverABM(s, dt=0.9*s.t_scale)

with NoInterrupt(ignore=True) as interrupted:
    while e.t < e.y.T and not interrupted:
        e.evolve(100)
        plt.clf()
        e.y.plot()
        display(plt.gcf())
        clear_output(wait=True)
```

# Solitons

+++

The GPE admits grey solitons with infinite period.  These solitons arise from the previous solutions in the limit where $m\rightarrow 1$ so that $\sn\rightarrow\tanh$ and $\cn\rightarrow\dn\rightarrow \sech$.

The soliton solution is:

$$
  \psi(x, t) = \sqrt{\bar{\rho}}e^{-\I\mu t/\hbar}\left[
  \I\frac{v}{c} 
  + \frac{u}{c}\tanh\left(\frac{x-vt}{l}\right)\right], \qquad
  \rho = \bar{\rho}\left[\frac{v^2}{c^2}
  + \frac{u^2}{c^2}\tanh^2\left(\frac{x-vt}{l}\right)\right]\\
  u^2+v^2 = c^2, \qquad \mu=g\bar{\rho} = mc^2,\qquad l=\hbar/m u.
$$

This can be expressed as

$$
  \psi(x, t) = \sqrt{\rho}e^{-\I\mu t/\hbar + \I\phi(x)}, \qquad
  \tan\phi(x) = \frac{v}{u\tanh\left(\frac{x-vt}{l}\right)}, \qquad
  \phi_{\mathrm{twist}} = \phi(\infty) - \phi(-\infty) = 2\tan^{-1}\frac{v}{u}
$$

+++

$$
  \rho(0) = \bar{\rho}\frac{v^2}{c^2}
$$

+++

or, in the same units $m=\hbar = g = 1$ as El and Hoefer:

$$
  \rho = \bar{\rho}\left[\frac{v^2}{c^2} + \left(1-\frac{v^2}{c^2}\right)\tanh^2[u(x-vt]\right]\\
  u= \sqrt{c^2-v^2}, \qquad \mu=\bar{\rho} = c^2
$$


To compare with (2.117):

$$
  \rho = \frac{(r_4-r_3-r_2+r_1)^2}{4} +(r_4-r_3)(r_2-r_1)\sn^2[\sqrt{(r_4-r_2)(r_3-r_1)}(x-Vt), m]
$$

we have $m=1$, $\sn = \tanh$ and

$$
  \frac{(r_2-r_1)(r_4-r_3)}{(r_4-r_2)(r_3-r_1)} = m = 1,\\
  \frac{(r_4-r_3-r_2+r_1)^2}{4} = \bar{\rho}\frac{v^2}{c^2} = v^2\\
  a = (r_4-r_3)(r_2-r_1) = \bar{\rho}\left(1-\frac{v^2}{c^2}\right) = \bar{\rho} - v^2\qquad
  \sqrt{(r_4-r_2)(r_3-r_1)} = u, \\
  V = \frac{r_1+r_2+r_3+r_4}{2}
$$

+++

Everything here is reasonable except the definition of the phase velocity $V$ which differs from the soliton velocity $v$.

+++

The background density is

$$
  \bar{\rho} = \frac{(r_4-r_3-r_2+r_1)^2}{4} + (r_4-r_3)(r_2-r_1)
             = \frac{(r_4-r_3+r_2-r_1)^2}{4}.
 $$

+++

To compare with (2.122):


, the stationary solution must have $\bar{\rho} = a_s$ or:

$$
  (r_4-r_1)^2 = 4(r_4-r_2)(r_2-r_1)
$$

If $r_1 = -r_4-2r_2$, then

$$
  4(r_4+r_2)^2 = 4(r_4-r_2)(r_4+3r_2)
$$

```{code-cell}

```

# Issues

+++

I was having some issues finding solutions so I considered stationary solutions:

$$
  \psi = \sqrt{\rho_0 + a\sn^2(kx,m)}.
$$

Symbolically solving for solutions to the GPE (in Maple) gives the following solutions:

\begin{align}
  a&=k^2m^2, &&& \mu &= \rho_0 - \frac{k^4m^2}{2\rho_0}, \\
  a&=0, &&& \mu &= \rho_0, & \psi &= \sqrt{\rho_0}.\tag{Constant}\\
  a&=\rho_0 m^2 , & k&=\sqrt{\rho_0}, & \mu &= \rho_0- \frac{\rho_0m^2}{2}, & \psi &= \sqrt{\rho_0}\sqrt{1+m^2\sn^2(\sqrt{\rho_0}x,m)}.\tag{Constant}\\  
  a&=\rho_0 , & k&=\frac{\sqrt{\rho_0}}{m}, & \mu &= \rho_0 - \frac{\rho_0}{2m^2}, & \psi &= \sqrt{\rho_0}\sqrt{1+\sn^2\left(\frac{\sqrt{\rho_0}}{m}x,m\right)}.\tag{Constant}\\  
\end{align}

(assuming $a,k\neq 0$ and real $k$) gives:

$$
  m^2 = \frac{a}{k^2}, \qquad
  \mu = \rho_0 - \frac{ak^2}{2\rho_0}, \qquad
  \rho = \rho_0[1 - \sn^2(kx,m)].
$$

```maple
restart;
X_:=JacobiSN(k*x,sqrt(m2));
psi:=sqrt(rho[0] + a*X_^2);
Rho:=psi^2:
res:=collect(numer(simplify(-diff(psi, x, x)/2/psi + Rho - mu)), X_):
mu:=expand(solve(coeff(res, X_, 0), mu));
m2:=solve(coeff(res, X_^6), m2);
solve([coeffs(res, X_)], [a,k^2]);
```

+++

The dark soliton solution has $m=1$, $a=k^2$

```{code-cell}
s.m
```

$$
  a = (r_4-r_3)(r_2-r_1)\\
  a = mk^2\\
  k^2 = (r_4-r_2)(r_3-r_1)\\
$$

```{code-cell}
s = State(rs=[-4.1, 1.0, 1.1, 2.])
x = s.xyz[0]
m = s.m_
a = s.a_
k = s.k_
rho_0 = s.rho_0_
print(m, a, m**2*k**2, rho_0, k)
k = np.sqrt(rho_0)
plt.plot(x, 1+sn(k*x,m)**2)
#s[...] = np.sqrt(rho_0*(1-sn(x,m)))
```

$$
  m^2 = \frac{a\alpha^2}{D^2}, \qquad
  (3m^2 - a) = 6\frac{a\alpha^2}{D^2}, \qquad
  -2a(m^2 + 1) = 6\frac{a\alpha^2}{D^2}, \qquad
  \mu = \frac{-D^2a}{2} + \alpha^2
$$

$$
  -2a(m^2 + 1) = (3m^2 - a) = 6m^2\\
  a = -3m^2\\
  m^2 + 1 = 1
$$

+++

To do this, we first consider boosting to a moving frame with velocity $v$ so that in this frame, the traveling wave solution is stationary and periodic.  We second 
Here we are looking for solutions that satisfy:

$$
  \psi(x + L) = e^{\I m v L / \hbar}\psi(x)
$$

```{code-cell}

```

# Minimization

+++

Here we consider the hypothesis that traveling waves can be found as minimum energy solutions in a box of period $L$ with twisted boundary conditions holding both the total particle number fixed and the value of the wavefunction at one point.

+++

We start with the soliton solution:

$$
  \rho = \bar{\rho}\left[\frac{v^2}{c^2} + \left(1-\frac{v^2}{c^2}\right)\tanh^2[u(x-vt]\right]\\
  u= \sqrt{c^2-v^2}, \qquad \mu=\bar{\rho} = c^2.
$$

At the core, the density is:

$$
  \bar{\rho}\frac{v^2}{c^2}
$$

+++

$$
  \psi(x, t) = \sqrt{n_0}\left[\frac{v}{c}\I + \sqrt{1-\frac{v^2}{c^2}}\tanh\left(\frac{x-vt}{l}\right)\right]e^{-\I\mu t/\hbar}\\
  \mu = gn = mc^2\\
  \psi_v(x_v) = \sqrt{n_0}e^{-\I mvx_v/\hbar}\left[\frac{v}{c}\I + \sqrt{1-\frac{v^2}{c^2}}\tanh\left(\frac{x_v}{l}\right)\right]e^{-\I(\mu + m v^2/2) t/\hbar}\\
  \mu = gn = mc^2  
$$

```{code-cell}
from gpe.imports import *
import gpe.Examples.traveling_waves;reload(gpe.Examples.traveling_waves)
from gpe.Examples.traveling_waves import (
    StateTravellingWave, MinimizeStateTravellingWave)
from gpe.exact_solutions import TravellingWaves

n1 = 1.0
n0 = 0.5
tw0 = TravellingWaves(n1=n1, n0=0.5)
n_avg = tw0.get_density().mean()
s = StateTravellingWave(Nx=tw0.basis.Nx, Lx=tw0.basis.Lx, 
                        psi_0=tw0[0], ind=0, 
                        n_avg=n_avg, v_p=tw0.v_x, twist=tw0.twist)
s.set_psi(tw0.get_psi())
#s[...] = ts0[...]
```

```{code-cell}
:init_cell: true

s.cooling_phase = 1
e = EvolverABM(s, dt=0.1*s.t_scale, normalize=False)
with NoInterrupt() as interrupted:
    while not interrupted:
        e.evolve(200)
        plt.clf()
        e.y.plot()
        display(plt.gcf())
        clear_output(wait=True)
```

```{code-cell}
plt.plot(tw0.xyz[0], s[...] - tw0[...])
```

```{code-cell}
plt.plot(tw0.get_Vext()- s.get_Vext())
```

```{code-cell}
dy0 = tw0.empty()
tw0.compute_dy_dt(dy=dy0, subtract_mu=False)
dy = s.empty()
s.compute_dy_dt(dy=dy, subtract_mu=False)
dy0[...] - dy[...]
```

```{code-cell}
s.plot()
```

```{code-cell}
dy = s.empty()
s.compute_dy_dt(dy=dy)
```

```{code-cell}
s.cooling_phase = 1j

e = EvolverABM(s, dt=0.5*s.t_scale, normalize=True)
with NoInterrupt() as interrupted:
    while not interrupted:
        e.evolve(100)
        plt.clf()
        e.y.plot()
        display(plt.gcf())
        clear_output(wait=True)
```

```{code-cell}

```

```{code-cell}
from scipy.optimize import brentq
n1 = 1.0
n_avg = 0.9
def f(n0):
    print(n0)
    tw = TravellingWaves(n1=n1, n0=n0)
    return tw.get_density().mean() - n_avg
brentq(f, 0, n1*0.999)
```

```{code-cell}
:init_cell: true

from gpe.imports import *
import gpe.Examples.traveling_waves;reload(gpe.Examples.traveling_waves)
from gpe.Examples.traveling_waves import StateTravellingWave, MinimizeStateTravellingWave
```

```{code-cell}
vs = np.linspace(-1,1,10)
Es = []
for v in vs:
    s0 = StateTravellingWave(Nx=128, L=10.0, v_p=v, twist=0.1)
    m = MinimizeStateTravellingWave(s0)
    s = m.minimize()
    Es.append(s.get_energy())
    plt.clf()
    s.plot()
    display(plt.gcf())
    clear_output(wait=True)
```

```{code-cell}
plt.plot(vs, Es)
```

```{code-cell}

```

Here is a stationary dark soliton.

```{code-cell}
%%time
dark_soliton = Soliton()

s0 = StateTravellingWave(psi_0=0, twist=np.pi, v_p=0)
s0[...] = dark_soliton.get_psi(s0.xyz[0])

m = MinimizeStateTravellingWave(s0)
s = m.minimize(E_tol=1e-12, psi_tol=1e-12)
s.plot()
plt.plot(x, abs(dark_soliton.get_psi(x))**2, '--', label='exact')
plt.legend(loc='best')
```

```{code-cell}
grey_soliton = Soliton(v=0.999)
plt.plot(np.angle(grey_soliton.get_psi(s0.xyz[0])))
```

```{code-cell}
v_c = 0.1
np.arctan2(v_c, np.sqrt(1-v_c**2))
```

```{code-cell}
v_c = grey_soliton.v/grey_soliton.c
theta = np.arctan2(v_c, np.sqrt(1-v_c**2))
s0 = StateTravellingWave(psi_0=0, twist=np.pi+2*theta, v_p=grey_soliton.v)
s0[...] = grey_soliton.get_psi(s0.xyz[0])
s0.psi_0 = abs(s0[...]).min()

m = MinimizeStateTravellingWave(s0)
s = m.minimize(E_tol=1e-12, psi_tol=1e-12)
s.plot()
plt.plot(x, abs(grey_soliton.get_psi(x))**2, '--', label='exact')
plt.legend(loc='best')
```

Here is a moving grey soliton.

```{code-cell}
plt.plot(np.angle(s._twist_phase_x)/np.pi)
```

```{code-cell}
%%time
rho_min = 0.01
v_p = np.sqrt(s.g*rho_min/s.m)
s0 = StateTravellingWave(psi_0=np.sqrt(rho_min), twist=np.pi, v_p=v_p)
s0[...] = s[...]
m = MinimizeStateTravellingWave(s0)
s = m.minimize(E_tol=1e-12, psi_tol=1e-12)
s.plot()
```

## Exact Solution

+++

$$
  \rho = \rho_{\min} + a \sn^2\left(\frac{x_v}{l};al^2\right), \\
  L = 2lK(al^2)\\
$$

```{code-cell}
from gpe.imports import *
#from mmfutils.math.special import ellipkinv
from scipy.special import ellipj, ellipk
import scipy.integrate
import scipy as sp
import gpe.bec

def sn(u, m):
    return ellipj(u, m)[0]

def K(m):
    return ellipk(m)


def get_QV(m=0.1, a=0.1, rho_min=1.0, N=256):
    rho_max = rho_min + a
    #m = a*l**2
    l = np.sqrt(m/a)
    L = 2*l*K(m)
    dx = L/N
    x = np.arange(N)*dx
    rho = rho_min + a*sn(x/l, m)**2
    C = np.sqrt(rho_min*rho_max*(1+rho_min*l**2))/l
    Vp = np.trapz(C/rho, x)/L
    Q = 2*np.pi/L
    return Q, Vp

rho_min = 1.0
a = 0.1
m = 0.63315
l = np.sqrt(m/a)
L = 2*l*K(m)
Q, Vp = get_QV(m=m, a=a, rho_min=rho_min)
```

```{code-cell}
x0 = m.pack(s0[...])
s0[...] = m.unpack(x0)
s0.plot()
```

```{code-cell}
s0 = State(Nx=64, L=L, mu=10.0, psi_0=np.sqrt(rho_min),
          v=Vp, m=1.0, hbar=1.0, g=1.0)
s0.plot()
plt.figure()
m = MinimizeState(s0, fix_N=True)
m.check()
s = m.minimize(psi_tol=1e-12, E_tol=1e-12)
s.plot()
```

```{code-cell}

```

```{code-cell}
s = State(Nx=128*8, L=120.0, mu=1.0, psi_0=0.5, v=v)
psi_s = s.exact_psi()
v = (np.angle(psi_s[-1]) - np.angle(psi_s[0])+np.pi)*s.hbar/s.m/(s.Lxyz[0] - s.Lxyz[0]/s.Nxyz[0])
s = State(Nx=128*8, L=120.0, mu=1.0, psi_0=0.5, v=v)
s[...] = s.exact_psi()
print(v)
```

```{code-cell}
plt.plot(s[...]/s.twist_phase)
```

```{code-cell}
v
```

```{code-cell}
s = s1
s.cooling_phase = 1.0
e = EvolverABM(s, dt=0.5*s.t_scale)

with NoInterrupt(ignore=True) as interrupted:
    while not interrupted:
        e.evolve(100)
        plt.clf()
        e.y.plot()
        display(plt.gcf())
        clear_output(wait=True)
```

```{code-cell}
s = State(Nx=128*4, L=20.0, mu=1.0, psi_0=0.1, v=0.1)
s[...] = 1.0
m = MinimizeState(s, fix_N=False)
s1 = m.minimize(use_scipy=True)
plt.plot(s1.xyz[0], s1.get_density() - abs(s1.exact_psi())**2)
#s1.plot()
abs(s1.get_density() - abs(s1.exact_psi())**2).max()
```

```{code-cell}
plt.plot(s1.xyz[0], np.angle(s1[...]))
plt.plot(s1.xyz[0], np.angle(s1.exact_psi()))
```

```{code-cell}
vs = np.linspace(0, 0.4, 10)
errs = []
for v in vs:
    s = State(Nx=128*4, L=20.0, mu=1.0, psi_0=0.2, v=v)
    m = MinimizeState(s, fix_N=False)
    s1 = m.minimize(use_scipy=True)
    errs.append(abs(s1.get_density() - abs(s1.exact_psi())**2).max())
```

```{code-cell}

```

# New

+++

$$
i \dot{\psi} = \left(E(\hat{k}) + gn \right)\psi, \qquad \psi(x,t) = \psi_v(x-vt)e^{i\phi(x,t)} 
= \psi_v(x-vt)e^{-i\mu t}
$$

$$
\left(E(\hat{k}) - v \cdot k + gn \right)\psi_v = 0
$$

```{code-cell}
from gpe.Examples import traveling_waves

class StateCustomDisp(traveling_waves.StateTwist_x):
    def __init__(self, Nx=128, L=10.0, mu=1.0, psi_0=1.0, ind=None,
                 v_p=0.0, twist=0.0, m=1.0, hbar=1.0, g=1.0):
        """
        Arguments
        ---------
        
        """
        self.v_p = v_p
        if ind is None:
            ind = Nx//2
        self.ind = ind
        self.psi_0 = psi_0
        self.p_v = p_v = m*v_p
        twist_x = twist  # + p_v * L / hbar
        traveling_waves.StateTwist_x.__init__(self, Nxyz=(Nx,), Lxyz=(L,),
                              m=m, hbar=hbar, g=g,
                              mu=mu, twist_x=twist_x)
        self[self.ind] = psi_0
    
    def init(self):
        traveling_waves.StateTwist_x.init(self)
        self._kx2 = (self.get_dispersion() * 2. * self.m / self.hbar**2)**2

    def get_Vext(self):
        return 0
        
    def get_dispersion(self):
        k = self._kx
        disp = (k * self.hbar)**2 / 2. / self.m - self.v_p * k  
        return disp
```

```{code-cell}
s.get_Vext()
```

```{code-cell}
s0 = StateCustomDisp(Nx=2**6, L=5, mu=1.0, psi_0=np.sqrt(0), v_p=0.1, twist=0)
m = MinimizeStateTravellingWave(s0)
s_init = m.minimize(E_tol=1e-12, psi_tol=1e-12)
```

```{code-cell}
def get_err(v_p, s_guess=[s_init]):
    s0 = StateCustomDisp(Nx=2**6, L=5, mu=1.0, psi_0=np.sqrt(0.1), v_p=v_p, twist=0)
    if s_guess[0] is not None:
        s0[...] = s_guess[0][...]
    m = MinimizeStateTravellingWave(s0)
    s = m.minimize(E_tol=1e-12, psi_tol=1e-12)
    s[s.ind] = s.psi_0
    s_guess[0] = s
    err = abs(s.compute_dy_dt(dy=s.copy())[...]).max()
    plt.clf()
    s.plot()
    plt.ylabel(v_p)
    plt.xlabel(err)
    display(plt.gcf())
    clear_output(wait=True)
    return err

v_ps = np.linspace(0,1,100)
errs = map(get_err, v_ps)
```

```{code-cell}
plt.plot(v_ps[1:], errs[1:])
```

```{code-cell}
v_p = 0.19
s0 = StateCustomDisp(Nx=2**6, L=10.0, mu=1.0, psi_0=np.sqrt(0.1), v_p=v_p, twist=0)
m = MinimizeStateTravellingWave(s0)
s = m.minimize(E_tol=1e-12, psi_tol=1e-12)
s.plot()
```

```{code-cell}
:comment_questions: false

sp.integrate.cumtrapz?
```

```{code-cell}
n = s.get_density()
L = s.basis.Lx
x = s.basis.xyz[0]
rho_min = n.min()
rho_max = n.max()

def get_rho_v_p(L=5.0, rho_min=0.1, rho_max=1.2, x=None):
    if x is None:
        x = np.linspace(-L/2., L/2.0, 1000)
    a = rho_max - rho_min
    m = get_m(a=a, L=L)
    l = np.sqrt(m/a)
    rho = rho_min + a*sn(x/l, m)**2
    C = np.sqrt(rho_min*rho_max*(1.+rho_min*l**2))/l
    plt.plot(x, rho)
    v_p = np.trapz(C/rho, x)/L
    u = v_p - C/rho
    theta = sp.integrate.cumtrapz(u, x, initial=0)
    return lambda x:rho_min + a*sn(x/l, m)**2, v_p, theta
```

```{code-cell}
plt.plot(x, theta)
```

```{code-cell}
plt.plot(x, s0[...].real)
plt.plot(x, s0[...].imag)
```

```{code-cell}
L = 3.0
rho_max = 1.2
rho_min = 0.5

rho, v_p, theta = get_rho_v_p(L=L, rho_min=rho_min, rho_max=rho_max)

s0 = StateCustomDisp(Nx=2**7, L=L, mu=rho_max, psi_0=np.sqrt(rho_min), v_p=v_p, twist=0)
x = s0.basis.xyz[0]
rho, v_p, theta = get_rho_v_p(L=L, rho_min=rho_min, rho_max=rho_max, x=x)
s0[...] = np.sqrt(rho(s0.basis.xyz[0]))*np.exp(1j*theta)
N_tot = s0.get_N()
Nx = s.basis.Nx
#s0[...] = np.sqrt((N_tot*Nx/L  - rho_min)/(Nx-1))*np.exp(1j*s.m*x*v_p/s.hbar)
#s0[s0.ind] = s0.psi_0
#print(N_tot, s0.get_N())


m = MinimizeStateTravellingWave(s0)
s = m.minimize(E_tol=1e-12, psi_tol=1e-12)
m.check()
s.plot()
#print s[s.ind]/s.psi_0 - 1
dy = s0.compute_dy_dt(s0.copy())
plt.plot(x, abs(dy[...]))
abs(dy[...]).max()
```

```{code-cell}
s[s.ind] = s.psi_0
    s_guess[0] = s
    err = abs(s.compute_dy_dt(dy=s.copy())[...]).max()
    plt.clf()
    s.plot()
    plt.ylabel(v_p)
    plt.xlabel(err)
    display(plt.gcf())
    clear_output(wait=True)
    return err

v_ps = np.linspac
```

```{code-cell}
s0 = StateCustomDisp(Nx=2**6, L=10.0, mu=1.0, psi_0=np.sqrt(1.), v_p=0.1, twist=0.0)
m = MinimizeStateTravellingWave(s0)
s = m.minimize(E_tol=1e-12, psi_tol=1e-6)
s.plot()
plt.figure()
dy = s.copy()
dy = s.compute_dy_dt(dy=dy)
dy.plot()
```

```{code-cell}
dy = s.copy()
dy = s.compute_dy_dt(dy=dy)
plt.semilogy(abs(dy[...]))
```

$$
  \psi(x, t) = \sqrt{\rho(x, t)}e^{\I\phi(x, t)}
  \psi'(x, t) = \left(
      \frac{\rho'(x,t)}{2\sqrt{\rho(x, t)}} 
      + \I\phi'(x,t)\sqrt{\rho(x, t)}
  \right)e^{\I\phi(x, t)}\\
  \psi''(x, t) = \left(
      \frac{\rho''(x,t)}{2\sqrt{\rho(x, t)}} 
      - \frac{[\rho'(x,t)]^2}{4\sqrt{\rho(x, t)}^3} 
      + \frac{2\I\phi'(x,t)\rho'(x,t)}{2\sqrt{\rho(x, t)}}
      + \I\phi''(x,t)\sqrt{\rho(x, t)}
      - [\phi'(x,t)]^2\sqrt{\rho(x, t)}
  \right)
  e^{\I\phi(x, t)}\\
  \psi''(x, t) = \left(
      \frac{\rho''(x,t)}{2\sqrt{\rho(x, t)}} 
      - \frac{[\rho'(x,t)]^2}{4\sqrt{\rho(x, t)}^3} 
      + \frac{2\I u(x,t)\rho'(x,t)}{2\sqrt{\rho(x, t)}}
      + \I u'(x,t)\sqrt{\rho(x, t)}
      - [u(x,t)]^2\sqrt{\rho(x, t)}
  \right)
  e^{\I\phi(x, t)}\\
$$  
$$
, \qquad
  u(x, t) = \phi'(x, t),\\
  m = al^2, \qquad
  \rho(x, t) = \rho_{\min} + a \sn^2\bigl(\frac{x_v}{l};m\bigr), \qquad
  x_v = x - Vt, \qquad
  u(x, t) = V - \frac{C}{\rho(x, t)}, \qquad
  \rho_{\max} = a + \rho_{\min}, \qquad
  C = \frac{\sqrt{\rho_{\min}\rho_{\max}(1+\rho_{\min}l^2)}}{l},\qquad
  L = 2lK(m), \qquad
  Q = \frac{2\pi}{L}.
$$

```{code-cell}

```

```{code-cell}

```